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Abstract. We investigate the properties of the stochastic Gross-Pitaevskii equationdescribing 
a condensate interacting with a stationary thermal cloud derived by Gardiner et al We 
find the appropriate Ehrenfest relations for the SGPE, including the effect of growth noise and 
projector terms arising from the energy cutoff. This is carried out in the high temperature 
regime appropriate for the SGPE of |2), which simplifies the action of the projectors. 

The validity condition for neglecting the projector terms in the Ehrenfest relations is found 
to be more stringent than the usual condition of validity of the truncated Wigner method 
or classical field method - which is that all modes are highly occupied. In addition it 
is required that the overlap of the nonlinear term with the lowest energy eigenstate of the 
non-condensate band is small. We show how to use the Ehrenfest relations along with the 
corrections generated by the projector to monitor dynamical artifacts arising from the cutoff. 

We also investigate the effect of using different bases to describe a harmonically trapped 
BEC at finite temperature by comparing the condensate fraction found using the plane wave 
and single particle bases. We show that the equilibrium properties are strongly dependent on 
the choice of basis. There is thus an optimal choice of plane wave basis for a given cut-off 
energy and we show that this basis gives the best reproduction of the single particle spectrum, 
the condensate fraction and the position and momentum densities. 



1. Introduction 

Since the experimental achievement of Bose-Einstein condensation |4], the theoretical tool 
of choice for describing condensates is the Gross-Pitaevskii equation (GPE) |5i, which has 
been found to describe a remarkably wide range of ultra-low temperature BEC physics. The 
extension of the GPE to finite temperatures in recent years |6| |8] [9] 1101 approaches the 
problem by treating the highly occupied modes of the condensate as a multimode 'classical 
field' . This has developed alongside various stochastic generalisations of the GPE based on 
the phase space methods of quantum optics II II [T21 1131 . The recent formulation of Gardiner 
and coworkers (T]|2| derives a stochastic Gross-Pitaevskii equation(SGPE) for the condensate 
band of a partially condensed Bose gas in contact with the non condensed thermal cloud. The 
treatment uses the truncated Wigner approximation (TWA) 11411151 which has been used with 
a certain degree of success to describe highly occupied optical fields for many years 1161 : 
however, even in quantum optics, where the occupation numbers of the modes of interest are 
enormous, this approach has its limitations. 

The validity of the TWA for partially Bose condensed gases has recently been 
investigated in some detail by Sinatra et al |3). The authors model a trapped finite temperature 
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Bose-Einstein condensate using the TWA and a plane wave basis; consequently, the method 
uses a momentum cut-off in the numerical representation of the quantum field. In the high 
temperature regime it was shown that the condition of validity for the TWA is that all modes 
must be highly occupied. The approach of Gardiner and Davis is to use a GPE that is 
projected into a basis that generates a strict energy cut-off for the harmonic trap, for which all 
modes are highly occupied in the low energy region. The weakly occupied modes above the 
cut-off are then traced out, generating noise terms in the stochastic Gross-Pitaevskii equation 
(SGPE). An important computational question then arises: What are the Ehrenfest equations 
for the SGPE, and under what conditions can we recover physical behaviour that resembles 
the familiar motion of the GPE? This will be addressed in the first part of this paper, where it 
is shown that when the spatial variation of the thermal cloud is neglected the GPE Ehrenfest 
equations can be extended to the SGPE. We derive exact expressions for the boundary terms 
in the Ehrenfest equations generated by the projector which we use to devise a method of 
assessing the influence of the projector on the dynamics - a method which we expect to be 
applicable under the conditions of the TWA. 

The second part of the paper considers the validity of using a plane wave basis to model 
a harmonically trapped condensate. Until very recently numerical methods based on plane 
wave representations have been favoured largely because of the speed advantage gained by 
using pseudo-spectral Fourier transform methods 1 17 |, and the conceptual ease of use that is 
inherent in the plane wave basis. However, there are a number of approximation issues that 
arise when using a plane wave basis to describe a trapped condensate at finite temperature 
which have not previously been investigated. We explore the link between the TWA or 
classical field approximation and the basis of representation in detail. To determine how 
best the plane wave basis may be used for a trapped BEC we compare the plane wave 
basis with an efficient spectral method based on numerical quadrature developed by Dion 
and Cances 1181 . and applied to the finite temperature GPE by Blakie and Davis 1191 . We 
show how to construct the optimal plane wave basis for a given cut-off, and verify that this 
basis gives the best agreement with the spectral basis when computing the condensate fraction 
for a partially condensed Bose gas. Variation of the basis generates significantly different 
condensate fractions, and increases the region of phase space that is not modeled accurately 
by the basis. In any representation this region of phase space must be significantly occupied 
due to the TWA validity condition, but it is minimized for the optimal plane wave basis. 

Although the spectral method is more difficult to implement and often more 
computationally expensive than plane wave approaches 1201 . there are clear advantages to 
using it to describe a trapped finite temperature condensate. The energy cut-off is well defined, 
and the high occupation condition can be imposed symmetrically in the phase space. There is 
another important feature that is basis dependent: the form of the Ehrenfest equation for the 
condensate band energy depends on the basis in which it is evaluated. Using the correct basis 
numerically thus becomes vital if one wants to obtain the correct dynamics for the condensate 
band of an open system using the TWA. 

2. Background 

It is well known that a solution of the GPE will satisfy the same Ehrenfest relations as 
solutions of the Schrodinger equation 1211 . We briefly reiterate these here to establish 
notation. 
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2.1. Ehrenfest relations for the Gross-Pitaevskii equation 

The GPE is the equation of motion for a field evolving according to the Hamiltonian functional 
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In this paper V(x, t) describes a general time dependent trapping potential, and u = 
Anh 2 a/m is the S-wave limit interaction strength |5'|. The quantities of interest are the energy 
density 
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where (A) = J dx ip*Ai[>, and the energy is simply (U) = H [ip, ip*]. We are working with 
the many particle wavefunction, so for simplicity of notation the condensate band occupation 
number will be written as (N) = J dx ip*ip. 

2.2. Properties of projectors 

Projectors that carry out the separation into upper and lower energy bands are defined by first 
separating the potential 

V(x,t) = V (x) + SV(x,t), (9) 

where the time invariant part is used to define the single particle Hamiltonian 

-h 2 v 2 



Ha 
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+ V (x), 



(10) 



and 6V(x, t) is arbitrary. The representation basis is provided by the eigenstates 

H n (x) = e n n (x). (11) 

We introduce orthogonal projection operators which are defined with respect to the single 
particle basis by their action on an arbitrary wavefunction V>(x) 



Mx) / rf 3 yC(y)^(y 

n<N c 



(12) 
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qv= E / d3 yC(y)^(y)> (13) 

n>N a 

and satisfy V + Q = 1 The index of summation n represents all the eigenvalues required 
to specify the complete set of single particle modes, and N c defines the cut-off energy. For 
convenience the cut-off is chosen so that N c as an integer; hereafter we will use the notation 

E -E ("> 

n<N c n 

for projected sums over the condensate band. 

We now restrict our attention to the condensate band itself so that t/j = Vijj for the 
wavefunctions of interest. The projected GPE (PGPE) corresponding to (0 can be written as 

^ = -j(l-C)W. (15) 
By writing the PGPE in terms of the projector orthogonal to V we can use the properties 

(Qtpy = Q*ip* (i6) 

y*d 3 x0*Q^= y^x^)*^ (i7) 

QVi) = Qip = QH ^ = (18) 

to extract the Ehrenfest relations for the condensate band wavefunction along with 
modifications that arise from the Q projector. We make use of these relations in section |4] 
We have also used the following notation for the complex conjugate projector 

Q>= C(x) f d 3 y &,(y)V(y). (19) 
The delta function restricted to the condensate band 

Mx,y) = 5>„(x)C(y) (20) 

n 

has the projection property 

J d 3 y 5 c (x,y)f(y) = P/(x) (21) 

for any function /. <5c(x, y) behaves like a true delta distribution for functions restricted to 
the condensate band: 

d 3 yfc(x,y)PV'(y) = ^(x), (22) 

which is equivalent to VP = V. Note that a straightforward application of dl7l shows that 
the PGPE is energy and number conserving. 



3. The stochastic Gross-Pitaevskii equation 

The SGPE formalism 1 2 1 separates the partially condensed system into a low energy subspace 
of modes, or condensate band, and its orthogonal complement , the union of which furnishes a 
complete basis. The non-condensate band is assumed thermalized, so that it may be described 
by Gaussian statistics and traced out. The non-condensate band thus plays the role of a thermal 
reservoir and acts as a damping mechanism for the condensate band, while the condensate 
band contains the condensate and its excitations, along with a low energy thermal component. 



Properties of the SGPE 



5 



While this description is somewhat more complicated that the PGPE 171 ISl l9l [TOI when 
treated in full, if we neglect the scattering term (which does affect the condensate band 
occupation) and take the limit of a broad thermal cloud, it can be reduced to a relatively 
simple equation of motion for the condensate band which is closely related to the PGPE. 
The complete derivation may be found in reference 0, but here we will briefly sketch the 
derivation, with a few minor changes of notation to make a transparent connection with the 
results to follow. 

To obtain the SGPE we proceed from the second quantised Hamiltonian for the system 
in the S-wave scattering limit 

H = J d 3 x*t(x) + V{x,t) S j #(x) + ~ J d 3 x * t (x)* t (x)^(x)*(x). (23) 

The field operator is split at the cut-off energy into *(x) = </>(x) + ^nc( x )> where the non- 
condensate field ?/>nc( x ) describes the high energy thermal modes. The commutator of the 
condensate band field operator is 

[^(x),0t(y)] =fc ( X)y ). (24) 

The thermal statistics of the non condensate field allow averages over many non-condensed 
field operators to be factorised and reduced to products of single particle Wigner functions 
.F^u, v). In this work we neglect the phase damping processes which lead to the 'scattering' 
terms in the master equation of 0. The growth/loss master equation for the reduced density 
matrix of the condensate band pc = TrNC (p) can be written in terms of the amplitudesf 

G<+)(U ' V ' £)= (2^ / d " Kl J ^ J d " K3 F K K i)^( u > K 2)[l + ^(u>K 3 )] 

X 8(wi +LU 2 -UJ 3 - e /fr) e -*(Ki+K 2 -K 3 ).v (25) 

and 

G ( -)(u,v,e) = e (e -^ /fcsT G (+) (u,v,e) (26) 

in the form| 

Pc|growth= / d 3 u J d 3 v[{G ( - ) (u,v,L c )^(u-v/2)}p Cj( /)t(u + v/2)] 

d 3 u J d 3 v [p c {G<-' (u, v, -L c )tf (u - v/2) } , <f>(u + v/2)] 

+ Jd 3 uJ d 3 v[{G (+) (u,v,-L c )0t( u -v/2)}p c ,0(u + v/2)] 

- Jd 3 u y'd 3 v[p c {G (+ »(u,v, J L c )0(u-v/2)},0t( u + v /2)] . (27) 
The condensate band operator is given in terms of the condensate band Hamiltonian 

h c = J d 3 x 0t (x) (-^ + y(x, t)j ^(x) + ^ y >x ^( x)0 t (x)0(x)0(x) (28) 

as 

L c ^(x) = [4>{x),H c ]. (29) 

f This corrects an extra minus sign in the defining equation (56) of (2j- 

| This corrects an error in equation (59) of wherein Lc appeared in place of — Lq in the second and third lines 
of the master equation equivalent to 1271 
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In principle a mean field or forward scattering term could also be included in He, and would 
alter the effective potential. In what follows we account for this possibility by including a 
general time dependent perturbing potential, into which such a term can be absorbed. 
The master equation for the growth is found by 

i) Expanding the exponential in the forward-backward relation \26\ to first order in 

(e - p)/k B T 

ii) Neglecting the condensate band energy during collisions with thermal atoms by 
approximating G <+) (x, v, e) « G <+) (x, v, 0) 

iii) Neglecting the condensate band momentum by making the approximation 0(u±v/2) w 
0(u) 

(i) is valid as long as the condensate band chemical potential is not significantly different 
from the non condensate chemical potential. This can be satisfied for a wide range of 
temperatures, and should be true for most physical situations of interest - the resulting SGPE 
is a valid description whenever the energy fluctuations between the two bands are small 
relative to k B T. A more detailed model of non condensate dynamics is required for this to be 
consistent in general since the dynamics of the non-condensate band should also be treated, 
but for many situations of interest the non condensate band can be described by a stationary 
distribution. Indeed, conditions (ii) and (iii) neglect energy and momentum exchange between 
the two bands, so the resulting equation is self consistent as long as the thermal cloud can 
be expected to remain stationary throughout the motion. Approximation (iii) can be made 
because G^ + ' (u, v, 0) is sharply peaked about v = 0. We further treat the growth amplitude 
G( + ' (u, v, 0) as spatially constant. This is the main simplifying approximation of the present 
work, and corresponds to the high temperature regime where the thermal cloud density is 
approximately constant over the condensate band. The growth parameter becomes 

7=7^^7^ / d 3 vG ,+, (0,v,0). (30) 



k B T k B T t 

Relaxing this assumption is possible in principle, but leads to Ehrenfest relations that have 
a complicated dependence on the shape of the thermal cloud and the cut-off energy of the 
projector. The form we use simplifies the projection as much as possible. 

The linearization required to obtain the SGPE from \21\ is an expansion in the operator 
j(Lc — fi)/k B T, which requires this to be small compared to the usual Gross-Pitaevskii 
evolution. The prefactor is usually of the order fry ~ 10~ 3 so this is easily satisfied in 
practice. In the high temperature regime y takes the simple form 

mk B Ta 3 

7 = z , (31) 

nu 

which we will use in this work. 

The SGPE derivation then follows the familiar path of quantum optics (see |2|): The 
master equation for the reduced density matrix of the condensate band is mapped onto an 
equation of motion for the multimode Wigner distribution. Derivatives higher than second 
order in the fields are neglected in order to derive a genuine Fokker-Planck equation with 
positive definite diffusion matrix. This is mapped to a stochastic differential equation for the 
condensate band field a(x, t) 1221 . The distinction between [;2i] and earlier work is a rigorous 
formulation in terms of projection operators that generate a consistent energy cut-off for the 
condensate band. 

We can now write the the SGPE as the Ito stochastic differential equation 

da(x, t) = -^L G pa(x, t)dt + jV(n - L G p)a(x, t)dt + dW G (x, t), (32) 
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where the noise is a vector Wiener process which satisfies 

dW G (x,t)dW G (x',t) = dW G (x,t)dW G (x',t) = (33) 
dW G (x,t)dW G (x' \t) = 2jk B T5 c (x,x')dt, (34) 

An important consequence of neglecting the scattering in 1 2 1 is that the growth noise is 
purely additive, so that (1321 is the same in Ito or Stratonovich form 1231 . 



3.1. The PGPE and Fudge equations 

Equation (I32> can be used to recover two useful equations which have already been used 
to successfully describe a number of interesting finite temperature effects in BEC including 
equilibrium properties and vortex lattice nucleation f7ll8ll9l lT0ll24ll25l . 
Putting 7 = recovers the PGPE (O 

ih ^&r~ = v { (~^§r + y(x ' t] + uHx ' * )|2 ) a(x ' t] } • (35) 

The projector is number conserving and acts to keep the wavefunction in a restricted energy 
subspace of the trap. The careful addition of this single feature to the GPE expands its scope 
to finite temperatures, where the field a(x, i) provides a non-perturbative description of both 
the condensate and a range of excitations up to the cut-off energy of the projector 1261 . 

A rigorous form of the phenomenological finite temperature GPE - or Fudge equation (2 
1241 - is found from d32t by dropping the noise and settting the projector to the identity, to give 



da(x f) / h 2 \7 2 \ 

ih — )-^- = (1-ihj) \-V(x,t) + u\a(x,t)\ 2 a(x,t) + ihyyua(x,t). (36) 

ot \ 2m J 

This is a semiclassical equation for a condensate in contact with a thermal cloud at chemical 
potential [i. Although many approximations have been made to derive this equation they are 
all well defined, and it must not be overlooked that this is therefore the logical equation to 
use when treating a damped condensate - rather than the phenomenological approaches to 
damping that have thus far been used 1271 . We also note that numerical integration of d36l > 
will lead to damping for any 7, rather than the small 7 limit required for the Fudge equation 
of 0. 

4. Continuity and Ehrenfest relations at finite temperature 

The projector can be dealt with by noting that, by putting V = 1 — Q in ( 1321 . the standard 
Ehrenfest and continuity results can be recovered from the terms multiplied by the identity, 
and the properties d!6i - (1181 can be used to find the influence of the Q projector. 

4.1. Continuity 

Averaging over the noise using Ito rules leads to the continuity equation 

+ V • j c (x) = ^Im(Q* {(SV(x,t) +u|a(x)| 2 )a*(x)}a(x)) u . 

+ 7 2Rc(Q* {[SV(x,t) +u\a(x)\ 2 ]a*(x)}a(x)) n . 

+ 7 2Re(/i|a(x)| 2 - a*(x)L GP a(x)) w (37) 
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where we have used and () w denotes the Wigner average. The density and current take 
the usual form in terms of the wavefunction a(x) and correspond to symmetrised operator 
averages for the condensate band 

n c (x) = i(^(x)0(x) + 0(x)^( x )) (38) 

jc(x) = ii(V0t(x)^(x) - t (x)V^(x) + V0(x)0t( x ) - 0(x)V0t(x)>, (39) 
4m 

where the angle brackets around an operator expression denote the trace over the density 
matrix of the condensate band. When there is no damping (7 = 0) the continuity equation 
for the resulting PGPE has an additional source term (the first term on the right hand side of 
equation \31V ) which redistributes the field. However, for any function /(x) and projected 
wavefunction a(x) 

j d 3 x Q*[/(x)]a(x) = J d 3 x /(x)Qa(x) = 0, (40) 

and consequently the source conserves atom number. The damping generates qualitatively 
different terms responsible for growth and fluctuations, as well as another number conserving 
source term. 



4.2. Ehrenfest relations 

The fluctuation terms have a generic form that can be written as a trace over the projected 
operator; the energy has to be treated slightly differently because of the nonlinear term, and 
may not be cast in the same form as the other averages. 

We have overloaded the () notation somewhat, since we wish to represent both spatial 
integration and stochastic averages to write down the Ehrenfest equations; in what follows we 
will adopt the convention 

(/),, = |d 3 x(a*(x)/(x)a(x)) wj (41) 

that is, we will suppress the coordinates of spatial integration except where they are required 
to specify the operators which are being integrated. To carry out the averaging, we use Ito 
rules to average over the noise, and the properties of the Q projector. In the equation of 
motion for (A), where A is one of the operators x, p or L, the Ito correction comes about 
from the last term in 

d(A) = J d 3 x a*(x)Ada(x) + da*(x)Aa(x) + da* (x)Ada(x), (42) 

Carrying out the spatial and stochastic averaging leads to the Ehrenfest relations for the 
projected stochastic Gross-Pitaevskii equation: 

d(N) w 



dt 

d(x)w = {p 
dt n 
d(p) w 



2 7 (/i - L GP ) W + 2 7 fc B TTr(n (43) 
h 2 7 Re(x( A i - L GP )) W + 2 1 k B TTr(Vx) + Q x , (44) 



dt 

d(L) w i 



(W) w + 2 7 Re(p( A/ - Lgp)) w + 2 1 k B TTr(Vp) + Q p , (45) 
LV) W + 7 2Rc(L( M - L GP )) lv + 2 7 fc B TTr(7>L) + Q L , (46) 



dt h 

+ 2 7 <T GP ( A1 - L GP )) W + 2 1 k B T(Tv[P(H + SV)] + u(S c )), 



dt \ dt 



(47) 
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where the left acting operator is 

a*(x)t GP = L GP a*(x), (48) 
and the boundary terms take the form 

Qa = hlm{F A ) w + 7 2Rc(F A ) w , (49) 

where 

F^a(x) = (5V(x,t) + u\a(x)\ 2 )Q[Aa(x)}. (50) 

These equations serve as useful consistency conditions for numerical simulations of the 
SGPE, and extend the Ehrenfest results to the projected formalism. The Ehrenfest behaviour 
familiar from the Gross-Pitaevskii equation is modified by growth terms and boundary 
corrections generated by the Q projector. 

It is worth noting that corresponding results for the Fudge equation d36l are obtained 
by putting T = in equations J43> - J47> . while retaining 7 and omitting the noise average. 
Although, from ( 13 11 we see that 7 is in fact proportional to T, the equations have been cast 
so that T only occurs explicitly in terms arising from the noise. It is clear from ( 149ft that 
when T = the boundary term caused by damping can still play a role in the Fudge equation 
dynamics. If we further put 7 = the results for the PGPE are recovered, which still have 
boundary corrections corresponding to the imaginary part of the source term in equation d37i . 

4.2.1. General remarks 
i) A typical trace term is, for example 



Tr(7>x) = Tr ^ |</»„)(0„|x = ^ f rf 3 x ? * (x)x0„(x) (51) 

n n J 

The assumption of a homogeneous non-condensate thus leads to state independent 
driving from the projected operators; the exception is the energy equation ( 147 1 , where 
the last term 

(6 C ) = J d 3 x a*(x)fc(x,x)a(x) (52) 

arises from the nonlinear interaction and cannot be cast as a trace. 
The term in i43i 

Tr (^) =l>«|0n) ( 53 ) 
n 

is just the number of modes in the condensate band. 

ii) Although all explicit projectors have been accounted for, the spatial integrals generate 
implicit projection since we are working with a projected stochastic wavefunction 
a = Va. 

iii) It is immediately apparent from ( I50t that the boundary terms vanish when either 
[A, H Q ] = 0, or 

Q* [(6V(x,t) + u|a(x)| 2 ) a*(x)] =0, (54) 

which is automatically true if the modes at the energy cut-off are weakly occupied. 
However, since the TWA has been used we have assumed that all modes in the condensate 
band are significantly occupied. These two conditions can be reconciled if the occupation 
of the modes at the cut-off is small relative to the other modes in the system but still large 
enough to render the third order derivatives in the Fokker-Planck equation unimportant. 
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iv) The advantage of neglecting the shape of the thermal cloud becomes clear when one 
considers the kind of terms that can show up in the spatially dependent calculation. 
A typical problematic term is proportional to QG(x)a(x), For a homogeneous cloud 
this term is zero, but when the spatial dependence becomes significant these terms, and 
additional ones that depend on AG(x) become important. 

v) Finally, we note the relationship between the number rate equation (I43> and the simple 
growth equation 1281 1291 

N Q = 2W + (N Q ) {(1 - e ^ c ^-"^ kBT )N Q + l} . (55) 

where Nq is the condensate occupation number, W + (Nq) is the growth rate into the 
condensate band, pc(No) is the condensate band chemical potential, and /i is the non- 
condensate band chemical potential. To the degree of approximation we are using for 
describing the non-condensate band we may linearize the exponential in (1551 and make 
use of 

W + (N a ) ee W + = yk B T (56) 
to give the linearized simple growth equation 

N = 2W + {lti~p c (N )]N /k B T + l}, (57) 
while, from d56i and J43i . we find the SGPE growth equation 

{N) w = 2W+ {(p - L G p) w Jk B T + M} , (58) 

where M = Tr V is the multiplicity of the condensate band. 

We note two main differences between the two growth equations d57t and J58l i: 

- The condensate chemical potential pc(No) of i51\ is replaced by a stochastic and 
spatial average of the Lq p operator over the condensate band. This is a significant 
generalization beyond the simple growth equation, and includes the influence of 
random initial conditions without making the random phase approximation required 
to obtain J55I . The great advantage of the SGPE approach is that it can describe 
both fluctuations and coherences within the condensate band, and may also be a 
good description for quasi condensate growth where phase fluctuations dominate 
the condensate band evolution. 

- A more subtle difference arises in the spontaneous terms in i57i (given by the 
additional +1 in braces). This therm is a consequence of treating the condensate 
as a single mode -that is, it is assumed to exist and be highly occupied |29|, and the 
rate equation describes those atoms in the condensate mode alone. 

The SGPE description is of the entire condensate band, and the condensate fraction 
must be extracted from the full field. The spontaneous term in d58l > (given by +M 
in braces) arises from the spontaneous growth into all of the single particle modes 
in the condensate band. 

4.2.2. Bose-Einstein condensate in a harmonic trap There is an important simplification of 
the Ehrenfest relations for a potential that is parity conserving. The operator traces that arise 
directly from the noise in d32l > are identically zero when the eigenstates conserve parity. This 
is easily seen for the parabolic trap, the case of most experimental interest. The driving terms 
in equations (1441 and ( 145 \ are proportional to Tr('Px) and Tr(Vp) respectively, consisting of 
sums of products of matrix elements of the form ((/> n \x\4> n ) , and (<f) n \p\<p n ) . These terms 
vanish when Vb(x) is parabolic because x and p couple harmonic oscillator eigenstates 
<t>n(x) — ► (j> n ±i( x )- Similar considerations give Tr('PL) = 0, as long as the thermal cloud is 
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stationary in the laboratory reference frame§. If the cut-off is chosen so that the Q projector 
terms are also small then we find the Ehrenfest relations for a trapped BEC 

d{N) w 



dt 
d{x) v 



2 7 ( M - L GP ) W + 2 1 k B T Tr(V) (59) 
2 7 Re(x(/x - L GP )) W (60) 



dt m 

= - (W) w + 2 7 Re(p(/i - L GP )) W (61) 
= - \(W) W + 2 7 Re(L( M - L GP )) W (62) 
' 2 1 (L GP {p - L GP )) W + 2 1 k B T(TT[P(H + SV)} + u(S c )) 



dt \ dt , 

(63) 

The fact that the finite temperature equations for (x) , (p) and (L) are the same as the equations 
found from J44l >— J46i by setting T = (which amounts to neglecting the noise) has interesting 
implications for condensate growth: the symmetry of the trap plays a crucial role in isolating 
the condensate band kinematic degrees of freedom from the direct influence of thermal noise. 
In a more general formulation fluctuations would still enter through the inhomogeneity of the 
thermal cloud, but these results are a good description in the high temperature regime. 

5. Projector terms for the harmonically trapped Bose gas 

One of the main aims of the classical field method is to deal with finite temperature BEC's, 
but the method is constrained by the requirement that all modes in the condensate band are 
highly occupied. It is clear that if there is significant occupation near the cut-off the dynamics 
can be radically altered, but it is not sufficient to simply monitor the occupation numbers. It 
is preferable to find a strict dynamical criterion that ensures the validity of the simulations. 
Rather than tackle this general problem at finite temperature, in this work we will simply 
show that at zero temperature it is possible to use the Ehrenfest relations to construct a 
reliable estimator of the error arising from the cut-off for a particular process: the Kohn mode 
oscillation. 

To consider a simple example of the phase space boundary effects (given by Qa in J43I— 
( I47l i) we consider a one dimensional model consisting of a harmonically trapped partially 
condensed Bose gas. We also take 7 = T = so that (N) = (E) = 0, andf 

dip 1 f d 2 ,\ , „ 
'it =2(~^ +X r' + ^ (64) 

d(x) 



dt 
dt 



(p) + Q x (65) 
- (x) + Q p , (66) 



where u is the dimensionless interaction strength. The details of how this relates to 
anisotropic trap geometry will not concern us here. We will merely consider modest effective 
nonlinearities of order 10 < Nil < 1000 to determine the validity of the criterion developed 
below. 

§ A treatment of rotating frame condensation will be carried out in 1301 

f We work in units of xq = {h/mnuj) 1 / 2 , to = a; -1 and fcrj = 1/xq for length and time and wave-vector 
respectively. 
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5.7. Harmonic oscillator basis 

We require some properties of the harmonic oscillator eigenstates. In terms of the Hermite 
polynomials H n (x) the eigenstates 

are coupled by the x and p operators to 

1 



c(f> n (x) = -j= (y/ruj> n -x{x) + \/ n + lcj) n+1 (x)) (68) 
■xj) n {x) = -j= (Vn+ l(f) n+ i(x) - y/n<pn-i(x)) , (69) 



so that the projector generates the terms 



Q xip(x) = \l c + 1 c Nc (t) cj) Nc+1 (x) (70) 



Qpip(x) = m/ ° + CAr c (t) 0at c+ i(x), (71) 



and the equations of motion are 

d{x 



(p) + Re B(N C + 1) (72) 
- (x) +JmB(N c + l), (73) 



dt 
dt 



B(n) = iuV2ncl_ 1 (t) / dx <fa(x) |^(a;)| 2 ^(a;). (74) 



where 



As expected the correction is essentially a boundary effect caused by the nonlinear term of the 
Gross-Pitaevskii equation: this is the only term in the Hamiltonian that can cause transitions 
between states in the condensate and non condensate bands. Since B(N C + 1) is proportional 
to c* N (t) weak occupation near the cut-off will naturally give Ehrenfest evolution. However, 
it is already apparent from the appearance of the nonlinear term that this condition alone is not 
sufficient to guarantee validity. We would like, therefore, to compute the relative error arising 
from the term B(N C + 1) in the above equations. Since the Ehrenfest equations involve 
easily computed averages, one could simply monitor the ratios |Re B(N C + l)/(p)\ and 
|Im B(N C + l)/(x)\. Unfortunately, these are unsuitable as estimators because (p) and (x) 
regularly pass through zero. However, we can exploit the rotational symmetry of the cut-off 
in phase space by using the complex variable (z) = (x) + i(p), for which the equation of 
motion is 

^ = - % ( z ) +B {N c + l). (75) 
dt 

We can ensure that the cut-off is not altering the dynamics appreciably by requiring 

B(N C + 1) 



< 1(T 4 . (76) 



To demonstrate this we use the following test, which we believe to be rather stringent: 
we take a ground state solution of the PGPE and give it a range of initial momenta in order 
to find the threshold where shape oscillations become apparent. For small kicks we observe 
Kohn mode oscillations and the projector has no effect on the dynamics. At a critical kick 
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(a) Density time series (b) PGPE ground state 



log(1+|W(x,k)|) log(1+|W(x,k)|) 




(c) Initial state shifted to fco = 8. (d) State after 1/4 trap period 



Figure 1. A simulation used to determine the threshold value of E z for a ground state 
wavefunction of the PGPE. (a) Distorted Kohn mode for motion near the semi-classical turning 
point of the cut-off mode (bold line), (b) Wigner function for the ground state solution of the 
fudge equation with nonlinear constant UN = 100, cut-off N c = 60. (c) Projected initial 
state after a shift in momentum space to wave vector fco = 8. (d) At 1/4 and 3/4 of a trap 
period the projector has its strongest influence. 



strength (determined by the cut-off energy), the edge of the condensate makes contact with 
the semi-classical turning point of the highest energy mode of the condensate band. The 
projector then comes into play, ensuring that the condensate wavefunction cannot make radial 
excursions in phase space that exceed the cut-off energy. This generates shape oscillations 
which become rather violent for large kick strength. A more revealing picture of the process 
is by transforming to phase space (which is detailed in the appendix), and the results for a 
high momentum simulation showing the distorted Kohn mode oscillations above threshold 
are presented in figure It is interesting to note that although the initial momentum kick 
has a negligible effect on the the shape of the condensate, the motion soon begins to fill 
the the available phase space - a kind of spurious thermalisation is caused by the interplay 
between the process of nonlinear mode mixing and the cut-off. The estimator we propose 
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+ Exact 




time [t J2n] 



Figure 2. The position and momentum averages (top) and their time derivatives (middle) show 
almost complete immunity to the projector, although small deviations are detectable at at one 
quarter and three quarters of a trap period. The cut-off error estimator is significantly larger 
than 10" 4 and has a well defined upper bound for a given initial energy. 



cut-off N c 


k 


max[iy/10- 4 


30 


1 


20 


40 


2.2 


10 


50 


2.5 


2 


60 


4 


2 


100 


7.3 


1 



Table 1. Threshold values for a range of cut-off energies. Other values are uN = 1 70, 
jU = 20. 



gives an unambiguous way of avoiding this predicament: it is clear from figure [2] that the 
lowest order moments themselves are rather insensitive to the cut-off effects, nevertheless, E z 
is an accurate measure of the artifact. 

We can now map the threshold in two different ways: we examine the dependence of the 
threshold value of E z on the cut-off at fixed nonlinearity in tabled and find the variation with 
nonlinearity for a given cut-off in table|2] The threshold value is independent of nonlinearity 
and cut-off energy over quite a wide range of either variable. It is important to note that, 
more generally, when any dynamical simulation is considered there will be a dominant set of 
moments that encapsulate the dynamics, and it is the effect of the projector on these quantities 
that must be considered. Higher order equations of the Ehrenfest type would then be used to 
extend the basic approach we have outlined here. 
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Nonlinearity uN 


ko 


max^J/10" 4 


59 


9 


2 


168 


6.7 


2 


476 


4.8 


8 


766 


3.3 


8 


1100 


2.2 


10 



Table 2. Threshold values for a range of nonlinear constants, with N c = 100. 



6. Comparison of projector effects for the plane wave and harmonic oscillator bases for 
the harmonic trap 

The classical field method requires that the majority of modes used in the simulation of a 
Bose gas are highly occupied. When imposing this condition consistently near equilibrium in 
a harmonic trap it becomes clear that a strict energy cut-off should be used. Such a cut-off is 
best implemented in the basis of harmonic oscillator eigenstates, since, for a well chosen basis, 
the full interacting Hamiltonian for the finite temperature system is approximately diagonal at 
the cut-off energy. Since the plane wave basis is often used for classical field simulations of 
trapped systems it is important to evaluate the validity of such a procedure against the more 
accurate procedure based on an exact energy cut-off. 

We can use the expressions for the effect of the projector on equations of motion for 
averages to get an idea of the artifacts that arise when using the plane wave basis to represent 
a trapped system at finite temperature. 



6.1. Plane wave basis 

The system defined by Vq(x, t) = SV(x, t) = 0, with periodic boundary conditions has been 
studied in great detail by Davis et al |9). We treat the one dimensional case. 



6.7.7. Plane wave basis without a trapping potential Writing 

k=-K 

where K is related to the spatial span according to K = 2ttN c / L. From this we have 

K 



which we can approximate by 



K 



J(k+A)x _ gi(ft— A)x 

2A : 



k=-K 

when the momentum grid spacing A is small. This leads to 

—i 



CKe i(K+A )x _ c _ Re -^K+A) x 



(77) 



(78) 



(79) 



(80) 



2AVL 

Since the momentum operator commutes with the Hamiltonian and the potential is absent, the 
equations for the operator averages are (TV) = (E) = (p) = and 

d(x) (p) 



dt 



+ Q X (K)+S«{L) 



(81) 
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where 



UK) 



-21m 



1 



dx |'^/'| 2 ' ! / , * 



2Av/L 



c K (t)e 



i(K+A)x 



-i(K+A)x 



sal) = ^ [{xi>)(j>r) - ripx^t'l, 

The boundary term S a arises from the finite span of the periodic basis, and will be important 
in situations where ip is non-zero at the box boundaries. This term does not arise in the use of 
the spectral basis, since the basis elements are defined over all space. 



(82) 



(83) 



6.7.2. Modeling a quadratic potential with a plane wave basis We have seen that there is a 
boundary term for the (x) equation when the plane wave basis is used for V(x,t) = 0. When 
the plane wave basis is used to model a system with a harmonic potential, the results follow 
from i50\ by putting 6V(x) = muj 2 x 2 /2, so that the entire potential becomes a variation with 
respect to H$. The continuity equation becomes 

+ V ■ j(x) = Him (Q* [(mwV/2 + u\^\ 2 )r] , (84) 
and the Ehrenfest relations become 

d{x) -^ + Qx (K) (85) 



dt m 

-muj' z (x) (86) 



dip) _ 



dt 

(87) 

with the boundary term given by 

Q X (K) = i 2Im( I dx [muj 2 x 2 /2 + u\^{x)\ 2 ) i/)*(x) (88) 

(89) 



2A-/L 



c K (t)e l{K+A > x - c_ K {t)e-^ K+A > 



We have neglected S a because the trapping potential will keep the wavefuction negligible 
at the grid edge for a well chosen basis. Comparing the plane wave basis Equations (1851 — 
J88i with the spectral basis results 1721 — 1741 (for SV = 0), we see that there are two 
notable differences: i) The boundary corrections occur in both the (x) and (p) equations 
for the spectral basis, and only in the (x) equation in the plane wave basis; and ii) there 
is a contribution from the potential in the boundary correction for (x) for the plane wave 
basis. This is particularly significant since this term can potentially assume large values when 
cx(t) ^ 0, even in the linear regime. 



6.2. The optimal plane wave representation 

There is a certain degree of freedom in choosing a plane wave basis for representing a 
harmonically trapped system. Here we show how to obtain the optimal plane wave basis 
that best captures the lowest harmonic oscillator states. We would expect this to be the best 
plane wave representation for modeling harmonically trapped systems. 

We consider a basis of N c plane wave states taken to extend over the spatial box of size 
x G [-L/2. L/2]. as defined in d77l . For fixed 7V C , the only free parameter in constructing the 
plane wave basis is L. Making L large is done at the expense of decreasing the momentum 
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width that can be represented on the grid, while conversely decreasing L limits the spatial 
extend of the system, but increases the momentum range. 

Here we give a simple argument for an optimal choice of L at fixed N c . In harmonic 
oscillator units the single particle Hamiltonian for the harmonic oscillator takes the form 



H=-k 2 
2 



1-2 

r 



(90) 



where bars are used to indicate dimensionless quantities (e.g. see Sec. VB1 of 1311 ). In these 
units the Hamiltonian and its eigenstates take the same form in coordinate and wave vector 
space. So the best grid choice will be when our numerical grids for (dimensionless) position 
and wave vector space are identical, i.e. when L = K. Returning to dimensioned units this 
optimal choice is 



L Q pt — 



2TlhN r 



(91) 



From this expression we directly obtain that the largest momentum the can be represented on 
the grid is i^cpt = TrN c /L opt (i.e. the limits of the sum in M1Y ). 

To quantify the sensitivity to non-optimal choices of L, in figure[3]we show the spectrum 
of energies found by diagonalizing the harmonic oscillator Hamiltonian in the plane wave 
basis for a range of L values. These results show that L opt is clearly the best choice, however 
even for L = L opt only about half the eigenstates are accurately obtained. 

6.3. Comparison of plane wave and harmonic oscillator phase space 




Figure 3. The numerical spectrum of the harmonic oscillator Hamiltonian. The solid lines are 
the plane wave results found by diagonalizing the harmonic oscillator Hamiltonian on a grid 
of N c = 16 points of width L. The dashed lines indicate the 16 lowest energies of the exact 
eigenspectrum, which corresponds to diagonalising the Hamiltonian in the spectral basis for 
N c = 16. 



The harmonic oscillator and plane wave bases differ somewhat in the regions of phase 
space they represent. This difference is reflected in the position quadrature grids associated 
with each basis, shown in figure|3a). It is apparent from this figure, that the spacing between 
quadrature points of the harmonic oscillator basis varies from being dense in the central region 
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to sparse at the edges. This enables the basis to capture large momentum states (i.e. fast 
spatial variations) at small displacements for the trap center, and smaller momentum states at 
large displacements. This suggests that for fixed N c , the harmonic oscillator basis captures a 
circular region of phase space, as shown schematically in figure^b). In contrast, the plane 
wave grid is equally spaced over the entire region [-L/2, L/2] (see figure^Ja)). This means 
the plane wave representation is equally well able to represent high momentum states at all 
displacements from the trap centre, suggesting that this basis captures a rectangular region of 
phase space (see figure|4lb)). 

Classically the motion of a harmonic oscillator corresponds to perfect circular trajectories 
in phase space, and we expect that a circular phase space projector forms the ideal energy cut- 
off. We note that the optimal plane wave representation corresponds to choosing L so that 
the maximum kinetic and potential energies associated with the edge values of the position 
and wave vector grids respectively are equal, i.e. ^mw 2 (L pt/2) 2 = h 2 (K opt /2) 2 /2m. For 
this case the phase space region bounded by the plane wave representations most closely 
matches the harmonic oscillator space (see the dashed line in figure |4jb)). In comparison, 
a non-optimal plane wave basis has an energy projector that restricts the kinetic energy and 
potential energy inconsistently, giving rise to a rectangular phase space boundary (e.g. see the 
dash-dot line in figure|4lb)). 

The high energy modes in the plane wave phase space exhibit anomalous dynamics 
arising from aliasing the region of phase space that is inconsistently represented. We illustrate 
this by examining the dynamics of a phase space point A in figure |4jb). This point evolves 
along the trajectory indicated by the arrow until it reaches the right position boundary (dashed 
line). It is then aliased to the left position boundary and continues to evolve through point 
B, before reaching the upper momentum boundary. It will then pass through point C, D and 
then return to A. The overall result is that states lying near the corners of the phase space 
region undergo a counter-clockwise evolution through phase space, in contrast to the normal 
clockwise evolution through phase space. In application to the SGPE formalism, we would 
expect that system disturbances with momentum and position characteristics lying in these 
corner regions will evolve in this anomalous manner. 

7. Comparison of plane wave and spectral representations for a thermalized ID gas 

In this section we compare the effect of basis on simulations of a harmonically trapped ID 
gas using the PGPE equation J35I . Our method follows the approach used by Davis et al. in 
reference (5): For each of the bases under consideration we evolve a randomized initial state 
of definite energy according to the PGPE. This evolution is expected to be ergodic, and by 
appropriately time averaging pure state expectations we are able to obtain ensemble averages. 
To compare the different bases we examine the equilibrium position and momentum density 
distributions, and the condensate fractions. 

In detail, the simulations we have conducted are for a dimensionless interaction strength 
of u — 200:z;o/cl> for bases with 40 modes (i.e. N c = 40). In Fig. |5]we present results for 
the density distributions found from evolving randomized initial states with a total energy of 
E = 14huj (Figs. |3a) and (b)), and for E = 2\ftw (Figs. |5Jc) and (d)). These two choices 
of energy correspond to a strongly condensed system (with a large condensate fraction) and 
a system close to the transition respectively (we will discuss condensate fraction later in this 
section). The three bases we compare are the spectral basis, and plane waves bases with 
L = L op t, L = 1.5L pt- 

The results in figure[5]confirm that the the optimal plane wave basis is in better agreement 
with the spectral bases than the non-optimal basis. A particular weakness of the plane wave 



Properties of the SGPE 



19 



( a ) Harm. osc. grid 

□ □□□□□□□□□□□□□□□□□□□□□□□□□□□□ □ 

Plane wave grid (Lopt) 

OOOCOCOCOOOOCOOC>CCOOOCXDOCCOCOO 

position [arb. units] 




-10 10 
position [arb. units] 



Figure 4. Phase space of the harmonic oscillator and plane wave representations, (a) The 
quadrature grids for the optimal plane wave basis (circles) and the harmonic oscillator basis 
(squares), (b) The approximate phase space captured by the optimal plane wave (dashed 
boundary), a non-optimal plane wave (dash-dot boundary), and harmonic oscillator (solid 
boundary) bases. The points A-D indicate the evolution of anomalous trajectories (see text). 




Legend: . . pw(1.5L ) 

1 opt' 

Harm. Osc. 



Figure 5. Results for the equilibrium position and momentum density profiles of a ID thermal 
Bose gas. Low Energy Case: Simulations for a total energy of E = lAtuo (a) equilibrium 
position density, (b) equilibrium momentum density. High Energy Case: Simulations for 
a total energy of E = 21Hlu (c) equilibrium position density, (d) equilibrium momentum 
density. Results are calculated by time-averaging classical field calculations carried out in 
the different bases under consideration: plane wave for L = L op t, L = 1.5Lopt, and the 
harmonic oscillator basis. Other simulation parameters: u = 200a.'o/i^, N c = 40. 



representations is that density distributions in the wings tend to decrease more slowly than 
the spectral basis. In the 1.5L opt case the density distribution even begins to increase near 
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Basis 


Cond. Frac. E = Uku 


Cond. Frac. E = 21huj 


H. Osc. 


0.370 


0.072 


PW (0.5L opt ) 


0.815 


0.238 


PW (0.8L opt ) 


0.392 


0.097 


PW(1.0L O pt) 


0.370 


0.078 


PW(1.2L opt ) 


0.360 


0.078 


PW(1.5L op t) 


0.444 


0.116 



Table 3. Condensate fraction for a ID thermal Bose gas obtained using various numerical 
bases. 

the boundary, as is apparent in figures |5J a) and (c). This behavior is most likely due to the 
inconsistent manner that the plane wave basis represents the highest energy states as discussed 
in the previous subsection. This will have serious implications for schemes that use the 
behaviour of the distribution wings to calculated temperature, and is likely to have effected 
the temperature calculations made in 1321 . 

Finally we examine how the condensate fraction is influenced by the choice of basis. 
We determine the condensate fraction using the Penrose-Onsager criterion 1331 1321 IT9l . To 
do this we calculate the 1 body density matrix by time-averaging the classical field, i.e. 
Pib(x,x') = (tjj* (x)ip(x'))timea,ve.- The condensate fraction is determined as the largest 
eigenvalue of the 1 body density matrix. The results are presented in table [5] for the cases 
considered in figure augmented by the results from a wider range of plane wave bases. 
These results show conclusively that non-optimal plane wave bases can have a dramatic 
influence on the physical properties of the system simulated. 

8. Conclusions 

We have shown how to derive exact Ehrenfest relations for the PGPE, an equation of 
motion which has become an important tool in the study of finite temperature Bose-Einstein 
condensates. The Ehrenfest relations show the interesting feature that for a harmonically 
trapped BEC in contact with a high temperature thermal cloud the operator averages for the 
kinematic degrees of freedom are immune to thermal noise on the average. 

For simulations of the PGPE and the SGPE in the truncated Wigner or classical field 
approximation it is required that all modes in the condensate band are significantly occupied 
0J|2JE1, however it is the relative occupation at the phase space boundary that determines 
the influence of the projector. Thus the cut-off can be chosen so that the projector does 
not generate spurious dynamics, even though the modes near the cut-off may have moderate 
occupation. However, this can be a delicate balance, and we have shown how to find easily 
computable dynamical tests of the PGPE that are not reliant on simply monitoring the mode 
occupation numbers. 

The projector in the SGPE generates boundary terms that arise because the Gross- 
Pitaevskii time evolution can evolve the wavefunction outside the condensate band. This 
kind of evolution can arise through either the nonlinear term or from an additional potential 
which is not part of the single particle Hamiltonian used to generate the representation basis. 
It therefore becomes important to choose the right basis, and we have shown that using the 
plane wave basis for a harmonically trapped BEC in thermal equilibrium can significantly 
alter the equilibrium condensate fraction unless the basis is chosen to optimally reproduce the 
single particle spectrum of the harmonic trap. 
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Appendix A. Transforming to phase space 

In order to examine the phase space behaviour of the condensate band wavefunction we 
require the connection between the mode representation ip(x) = ^„ c„ </>„ (x) and the Wigner 
function 

W(x, k ) = ^J d y e tky r(x + y/2)iP(x - y/2) (A.l) 

for the harmonic oscillator basis. 

We insert the mode decomposition into JA.U to get 

W(x, k)=J2 < c ™ W nm (x, k), (A.2) 

n,m 

where the modes 4> n (x) are the orthonormal eigenstates of the harmonic trap and 

W nm (x,k) = / dy e iky 4> n (x + y/2)4> m {x - y/2). (A.3) 



2tt _ 

For brevity we will use the notation 

W«(x,k) = Wn,n+ q {x,k). (A.4) 

We make use of the contour representation of the Hermite polynomials 1341 to write the 
modes as 

e -x 2 /2 r fa e -t 2 +2tx 

4> n (x) = ; <f> n (A.5) 

and (IA.4> becomes 

e~ x2 InHn + qV f ds e ~ s2 + 2sx 

r u p -t 2 +2tx 

D ~y 2 /4+y(ik+t- s ) 



, dy e~ y /i+v^+ts) (A.6) 
where the contours enclose the origin. Carrying out the y and t integrals gives 

-x 2 -k 2 / |9 9 1 f ^„ 2s(x+ik) 

Wil,.k) = L_^_JlH__^^L_ (a; _ ih _ sr+ ,, (A . 7) 

which, after the change of variables s = z(x — ik)/(z — 1), becomes 



W«(X, k) = ^4 Jt^tO* - i 4^T * ■ (A.8) 



7T 



(?i + <?)! y 2ttz / z n+1 {l-z) 



We can now use the generating function for the associated Laguerre polynomials 1341 

z) °° 

ZI = '£z n L q n (x) (A.9) 



e -xz/(l-z) 



(1 - 
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to find 



W*{x,k) = ^J-J-^—e-* ~ k (x - ikfLl{2{x 2 + k 2 )). (A.10) 

When g = 0we recover the Wigner function for a number state which is a well known result 
in quantum optics 1221 . 

Using the symmetry W n +q,n( x i k) = W%(x, k)*, the transformation to phase space is 

{N N-q ~\ N 

£ £ c* n c n+q WZ(x,k) - Yl \cn\ 2 W%{x,k). (A. 11) 

q=Q n=0 ) n=Q 
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